pckage <- c("polycor", "ggplot2", "gridExtra", "texreg", "lavaan",
            "ordinal", "tidyverse", "AER", "survey",
            "xtable", "nnet", "sandwich", "MASS", "lmtest", "semTable")
lapply(pckage, require, character.only=T)

#Set WD
setwd("C:/Users/Matias/Dropbox/Proyectos/Corruption and Political Knoweldge/Replication File")

#Open data with both waves merged in wide format
load("merged_ola_1y2.Rdata")
ls()
dim(ola)

##DROP ALL RESPONDENTS NOT INTERVIEWED IN WAVE 1
ola <- subset(ola, attrition==0)
dim(ola)


###################################################
############ SCALE CREATION #######################
###################################################

#------ Political Knowledge Index--------####

#Dimensionality of PK in Wave 1
pk_o1 = data.frame(with(ola, cbind(pk1_o1,pk2_o1,pk3_o1,pk4_o1,pk5_o1)))
pk_o1F <- pk_o1
for(i in 1:dim(pk_o1F)[2]) pk_o1F[,i] = as.factor(pk_o1F[,i])
factanal(cov=hetcor(pk_o1F, ML=F)$cor, factors=1)

#Ordinal alpha PK Wave 1
pk_o1_cor <- hetcor(pk_o1F, ML=F)$cor
mean_rho_pk_o1 <- mean(as.vector(pk_o1_cor)[as.vector(lower.tri(pk_o1_cor))==TRUE])
#Ordinal Alpha
(dim(pk_o1_cor)[1] * mean_rho_pk_o1) / (1 + mean_rho_pk_o1*(dim(pk_o1_cor)[1]-1))

#Dimensionality of PK in Wave 2
pk_o2 = data.frame(with(ola, cbind(pk1_o2, pk2_o2, pk3_o2, pk4_o2, pk5_o2)))
pk_o2F <- pk_o2
for(i in 1:dim(pk_o2F)[2]) pk_o2F[,i] = as.factor(pk_o2F[,i])
factanal(cov=hetcor(pk_o2F)$cor, factors=1)

pk_o2_cor <- hetcor(pk_o2F)$cor
mean_rho_pk_o2 <- mean(as.vector(pk_o2_cor)[as.vector(lower.tri(pk_o2_cor))==TRUE])
#Ordinal Alpha
(dim(pk_o2_cor)[1] * mean_rho_pk_o2) / (1 + mean_rho_pk_o2*(dim(pk_o2_cor)[1]-1))

#Political Knowledge Index 
ola$f1 <- apply(pk_o1, 1, sum, na.rm=T)
ola$f2 <- apply(pk_o2, 1, sum, na.rm=T)


#------Create Poltiical trust Index--------####

trstpol_o1 <- with(ola, cbind(trParlamen_o1, trGob_o1, trPartidos_o1, trJueces_o1, 
                              trFiscal_o1))
trstpol_o2 <- with(ola, cbind(trParlamen_o2, trGob_o2, trPartidos_o2, trJueces_o2, 
                              trFiscal_o2))
ola$pol_trst_o1 <- apply(with(ola, cbind(trParlamen_o1, trGob_o1, trPartidos_o1,
                                         trJueces_o1,trFiscal_o1)), 1, mean, na.rm=T)
ola$pol_trst_o2 <- apply(with(ola, cbind(trParlamen_o2, trGob_o2, trPartidos_o2,
                                         trJueces_o2, trFiscal_o2)), 1, mean, na.rm=T)
#Ordinal alpha Pol Trsut Wave 1 & 2
pt_o1_cor <- hetcor(trstpol_o1, ML=F)$cor
mean_rho_pt_o1 <- mean(as.vector(pt_o1_cor)[as.vector(lower.tri(pt_o1_cor))==TRUE])
#Ordinal Alpha
(dim(pt_o1_cor)[1] * mean_rho_pt_o1) / (1 + mean_rho_pt_o1*(dim(pt_o1_cor)[1]-1))

pt_o2_cor <- hetcor(trstpol_o2, ML=F)$cor
mean_rho_pt_o2 <- mean(as.vector(pt_o2_cor)[as.vector(lower.tri(pt_o2_cor))==TRUE])
#Ordinal Alpha
(dim(pt_o2_cor)[1] * mean_rho_pt_o2) / (1 + mean_rho_pt_o2*(dim(pt_o2_cor)[1]-1))


###########################################################
######------- FINAL MODELS WITH ROBUST SE' -----###########
###########################################################

#### TABLE 1 #####

# Linear Models - Cross-Lagged via OLS
m1a <- lm(corrup_comp_o2 ~ corrup_comp_o1 + f1 + lrscaleR_o1 + SEXO_PS.x + EDAD_PS.x +
            educa_o1F, data=ola)
m1b <- lm(f2 ~ corrup_comp_o1 + f1 + lrscaleR_o1 + SEXO_PS.x + EDAD_PS.x +
            educa_o1F, data=ola, subset=!is.na(corrup_comp_o2))

se_m1a <- sqrt(diag(vcovHC(m1a, type="HC0")))
se_m1b <- sqrt(diag(vcovHC(m1b, type="HC0")))
pval_m1a <- (1-pnorm(abs(m1a$coef / se_m1a)))*2
pval_m1b <- (1-pnorm(abs(m1b$coef / se_m1b)))*2

htmlreg(file="Table1.doc", list(m1b, m1a),
        override.se = list(se_m1b, se_m1a),
        override.pvalues = list(pval_m1b, pval_m1a), single.row = F,
        stars = c(0.01, 0.05, 0.1),
        caption="Table 1: Cross-Lagged Panel Model for Political Knowledge and Perceived Corruption",
        custom.coef.map=list("(Intercept)"="Intercept", "corrup_comp_o1"="Percieved Corruption W1",
                             "f1"="Political Knowledge W1",
                             "lrscaleR_o1Center"="Center W1",
                             "lrscaleR_o1Right"="Right-wing W1",
                             "lrscaleR_o1None"="No Ideology W1",
                             "SEXO_PS.xMujer"="Female W1",
                             "EDAD_PS.x"="Age W1",
                             "educa_o1FSecondary"="Secondary Ed. W1",
                             "educa_o1FTertiary Technical"="Tertiary Technical Ed. W1",
                             "educa_o1FTertiary Universitary"="Tertiary Universitary Ed. W1"),
        custom.model.names=c("Pol. Know. W2", "Per. Corrup. W2"))

screenreg(list(m1b, m1a), override.se = list(se_m1b, se_m1a),
          override.pvalues = list(pval_m1b, pval_m1a), single.row = F,
          stars = c(0.01, 0.05, 0.1))
        
#### TABLE 2: ROBUSTNESS & INTERACTION #####

m2a <- lm(f2 ~ corrup_comp_o1 + corrup_comp_o2 + f1 + lrscaleR_o1 + SEXO_PS.x + EDAD_PS.x +
            educa_o1F, data=ola)
m2c <- glm(f2 ~ corrup_comp_o1 + f1 + lrscaleR_o1 + SEXO_PS.x + EDAD_PS.x +
             educa_o1F, data=ola, subset=!is.na(corrup_comp_o2), family = quasipoisson)
m2d <- glm(f2 ~ corrup_comp_o1 + corrup_comp_o2 + f1 + lrscaleR_o1 + SEXO_PS.x + EDAD_PS.x +
             educa_o1F, data=ola, subset=!is.na(corrup_comp_o2), family = quasipoisson)
m3a <- lm(f2 ~ corrup_comp_o1 * lrscaleR_o1 + f1 + SEXO_PS.x + EDAD_PS.x +
     educa_o1F, data=ola, subset=!is.na(corrup_comp_o2))

se_m2a <- sqrt(diag(vcovHC(m2a, type="HC0")))
se_m2c <- sqrt(diag(vcov(m2c)))
se_m2d <- sqrt(diag(vcov(m2d)))
se_m3a <- sqrt(diag(vcovHC(m3a, type="HC0")))

pval_m2a <- (1-pnorm(abs(m2a$coef / se_m2a)))*2
pval_m2c <- (1-pnorm(abs(m2c$coef / se_m2c)))*2
pval_m2d <- (1-pnorm(abs(m2d$coef / se_m2d)))*2
pval_m3a <- (1-pnorm(abs(m3a$coef / se_m3a)))*2

htmlreg(file="Table2.doc", list(m2a, m2c, m2d, m3a), 
        override.se = list(se_m2a, se_m2c, se_m2d, se_m3a),
        override.pvalues = list(pval_m2a, pval_m2c, pval_m2d, pval_m3a),
        stars = c(0.01, 0.05, 0.1), single.row = F,
        caption="Table 2: Linear and Quasi-Poisson Panel Model for Political Knowledge",
        custom.coef.map=list("(Intercept)"="Intercept", "corrup_comp_o1"="Percieved Corruption W1",
                             "corrup_comp_o2"="Percieved Corruption W2",
                             "f1"="Political Knowledge W1",
                             "lrscaleR_o1Center"="Center W1", 
                             "lrscaleR_o1Right"="Right-wing W1",
                             "lrscaleR_o1None"="No Ideology W1", 
                             "SEXO_PS.xMujer"="Female W1",
                             "EDAD_PS.x"="Age W1",
                             "educa_o1FSecondary"="Secondary Ed. W1", 
                             "educa_o1FTertiary Technical"="Tertiary Technical Ed. W1", 
                             "educa_o1FTertiary Universitary"="Tertiary Universitary Ed. W1",
                             "corrup_comp_o1:lrscaleR_o1Center"= "PC W1 * Center W1",
                             "corrup_comp_o1:lrscaleR_o1Right"= "PC W1 * Right W1",
                             "corrup_comp_o1:lrscaleR_o1None"= "PC W1 * No Ideology W1"),
        custom.model.names=c("Model 3: Linear", "Model 4: Quasi-Poisson", 
                             "Model 5: Quasi-Poisson", "Model 6: Linear"))


screenreg(list(m2a, m2c, m2d, m3a), 
          override.se = list(se_m2a, se_m2c, se_m2d, se_m3a),
          override.pvalues = list(pval_m2a, pval_m2c, pval_m2d, pval_m3a))

##-----Complementary Analysis in Text----##

#Difference in predicted values of poisson models
nd <- with(ola, data.frame(corrup_comp_o1=c(4, 5), f1=mean(f1), lrscaleR_o1="Left", 
                           SEXO_PS.x="Mujer",
                           EDAD_PS.x=mean(EDAD_PS.x), educa_o1F="Secondary"))
predict(m2c, nd, type="response")


## Significance of effects of perceived corruption for each ideological group
#Centrist
beta <- m3a$coef
cov <-  vcovHC(m3a, type="HC0")
moderator <- c("lrscaleR_o1Center", "lrscaleR_o1Right", "lrscaleR_o1None")

#Effects
eff <- sapply(moderator, function(x) beta["corrup_comp_o1"] +
                beta[paste0("corrup_comp_o1:", x)])
#SE Effects
se_eff <- sqrt(sapply(moderator, function(x) {
  cov["corrup_comp_o1", "corrup_comp_o1"] +
    cov[paste0("corrup_comp_o1:", x), paste0("corrup_comp_o1:", x)] + 
    2*cov["corrup_comp_o1", paste0("corrup_comp_o1:", x)]
}))

#Pvalue
round((1-pnorm(abs(eff / se_eff)))*2,3)